% ! problem10.m --> based on <<MATLAB codes for Finite Element Analysis>>


clear;
clc;
format compact;

p=struct();

% define model paraeters
p.E=210e6;
p.L=10;
t=p.L/1000;
p.Iz=(1*t^3)/12;
p.P=-1000;
% define elems of bornoulli beam
p.elem_num=20;  
% dof of each nornoulli beam nodes
p.bornoulli_node_dof_num=2;
p.spring_k=100;

% -------------------------------

% node's coordinate
p.nodes=linspace(0,p.L,p.elem_num+1)';
p.elems=zeros(p.elem_num,2);
for i=1:p.elem_num
    p.elems(i,1)=i;
    p.elems(i,2)=i+1;
end
p.node_num=size(p.nodes,1);
p.node_Coord_x=p.nodes(:,1);

% global degree of freedom number: 
%   ! bornoulli_beam node's dof + 1 spring node's dof
p.global_dof_num=2*p.node_num;

% define global Nodal displacement  colum vector
p.displacements=zeros(p.global_dof_num+1,1);
p.node_forces=zeros(p.global_dof_num+1,1);

% define boundary conditions
%   ! BC: clamped at x=0
p.fix_dof=[1,2,p.global_dof_num+1]';

% initial global stiffness matrix
p.global_stiffness_matrix=zeros(p.global_dof_num+1);

% compute all ElemStiffnessMatrix and assembly stiff matrix
beam_dof_num=p.bornoulli_node_dof_num*size(p.elems,2);
p.elemStiffs=zeros(beam_dof_num,beam_dof_num,p.elem_num);

for i=1:p.elem_num
    connectivity=p.elems(i,:);
    % elem's all dof : a 3d truss node has 2 node,each node have 2 dof,
    % i-th node --> w_(2*i-1) / dw_(2*i-2) /  dof(displacement of node) 
    elem_dof=[connectivity(1)*2-1 connectivity(1)*2  ...
        connectivity(2)*2-1 connectivity(2)*2 ];

    node1_x=p.node_Coord_x(connectivity(1));
    node2_x=p.node_Coord_x(connectivity(2)); 
    
    elem_length=node2_x-node1_x;
    a=0.5*elem_length;
    k_e=(p.E*p.Iz/(a^3))*0.5*[3 3*a -3 3*a;
                        3*a 4*a*a -3*a 2*a*a;
                        -3 -3*a 3 -3*a;
                        3*a 2*a*a -3*a 4*a*a];
    p.elemStiffs(:,:,i)=k_e;
    p.node_forces(elem_dof)=p.node_forces(elem_dof)+((p.P*a/3)*[3 a 3 -a])';

    % assemble stiffness matrix
    p.global_stiffness_matrix(elem_dof,elem_dof)=...
        p.global_stiffness_matrix(elem_dof,elem_dof)+k_e;
end

% add spring nodes into global_stiffness_matrix
p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])=...
    p.global_stiffness_matrix([p.global_dof_num-1,p.global_dof_num+1],[p.global_dof_num-1,p.global_dof_num+1])+p.spring_k*[1,-1;-1,1];

% apply boundary condition
p.node_forces(p.fix_dof)=0;

% solution KU=F
p=solutionStruct(p);

% ! display result
disp('Bornoulli beam result !')

disp('beam挠度:')
omega_index=(1:2:p.global_dof_num)';
disp([omega_index,p.displacements(omega_index)])

disp('beam转角:')
theta_index=(2:2:p.global_dof_num)';
disp([theta_index,p.displacements(theta_index)])

disp('beam支反力:')
disp([p.fix_dof,p.node_forces(p.fix_dof)])

plot(p.node_Coord_x,p.displacements(omega_index),'-')
title('Problem10.m')